#usr/bin/env python
#
###################################################################
#
#module name: r3.strata for GRASS 7.0
#Author: Tim Bailey, Humboldt State University
#       
#Purpose: To populate voxel maps using data from profile and horizon data
#
#Copyright: 2013
#
#
#This program is free software under the GNU general Public License (>=v2)
#Read the file COPYING that comes with GRASS for details
#
####################################################################

#%module: r3.strata
#%description
#%keywords:rast3
#%keywords:voxel
#%keywords:horizon
#%keywords:profile
#%end
#%option
#%  key: rast3in
#%  type: string
#%  gisprompt: rast3 input
#%  description: input rast3 voxel map of lower boundaries of horizons
#%          This map can be generated from a 3d vector map with v.to.rast3
#%          v.to.3d is also potentially useful
#%  required: yes
#%end
#%option
#%  key: rast3out
#%  type: string
#%  description: output rast 3 voxel map
#%  required: yes
#%end
#%option
#%  key: doixy
#%  type: int
#%  description: distance of influence in the xy directions, number of cells that an identity propagates
#%  required: yes
#%end
#%option
#%  key: doiz
#%  type: int
#%  description: distance of influence in the z direction
#%  required: yes
#%end


import sys
import os
import grass.script as grass
#import atexit      not imlemented yet


def main():
    rast3_in= options['input']                  #input r3 map
    rast3_inout= options['output']              #output r3 map
    #doixy= int(options[doixy])                  #distance of influence in xy, number of cells
    #doiz= int(options[doiz])                    #distance of influence in z, number of cells
    #doixy= 20
    #doiz= 10

    def profile_builder():
#           n= global(doiz)   getting the doiz value into the iteration range, not implemented correctly
        for s in xrange(1, 10):
            b = s + 1
            m = 'z_iteration' + str(b)
            print m
            grass.mapcalc("profile_out = if(profile_in > 0, profile_in, if(profile_in[0,0,-1] > 0, profile_in[0,0,-1], 0))")
#profile cell values are unchanged or
#propagate upwards into null cells, one cell per iteration.
#actually in this version nulls turned to zero.
            grass.run.command('g.rename, input=profileout, output=profilein, --overwrite')
#reset the profile voxel map for mapcalc

    def horizon_builder():
#       for s in range(doixy)
        for s in xrange(1, 20):
            b= s + 1
            m= 'xy_iteration' + str(b)
            print m 
            vx_in = 'profile_in'#final iteration of profile builder is input for horizon builder routine
            grass.mapcalc("vx_out = if(vx_in > 0, vx_in,if(vx_in[1,0,0] > 0, vx_in[1,0,0], if(vx_in[0,1,0] > 0, vx_in[0,1,0],if(vx_in[-1,0,0] > 0, vx_in[-1,0,0], if(vx_in[0,-1,0] > 0, vx_in[0,-1,0])))))")
#Cell values propogate out in four directions from each profile voxel
            grass.run.command('g.rename, input=vx_out, output=vx_in --overwrite')#reset the voxel map name for mapcalc

    grass.run.command('g.rename, input=vx_in, output= output --overwrite')

if __name__ == "__main__": 
    options, flags = grass.parser() 
    main()

